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ABSTRACT 

Most numerical investigations on the role of magnetic fields in turbulent 
molecular clouds (MCs) are based on ideal magneto- hydrodynamics (MHD). 
However, MCs are weakly ionized, so that the time scale required for the mag- 
netic field to diffuse through the neutral component of the plasma by ambipolar 
diffusion (AD) can be comparable to the dynamical time scale. We have per- 
formed a series of 256 3 and 512 3 simulations on supersonic but sub-Alfvenic 
turbulent systems with AD using the Heavy-Ion Approximation developed in 
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Li et al.l (120061 ). Our calculations are based on the assumption that the number 
of ions is conserved, but we show that these results approximately apply to the 
case of time-dependent ionization in molecular clouds as well. Convergence stud- 
ies allow us to determine the optimal value of the ionization mass fraction when 
using the heavy-ion approximation for low Mach number, sub-Alfvenic turbulent 
systems. We find that ambipolar diffusion steepens the velocity and magnetic 
power spectra compared to the ideal MHD case. Changes in the density PDF, 
total magnetic energy, and ionization fraction are determined as a function of 
the AD Reynolds number. The power spectra for the neutral gas properties of 
a strongly magnetized medium with a low AD Reynolds number are similar to 
those for a weakly magnetized medium; in particular, the power spectrum of the 
neutral velocity is close to that for Burgers turbulence. 

Subject headings: MHD — turbulence — ISM: magnetic fields — ISM: kinematics 
and dynamics — methods: numerical 



Introduction 



Both supersonic turbulence and magnetic fields are widely observed in molecular clouds 
(MCs). M Cs have broad line width s, ranging from a few to more than 10 times the sound 
speed, c s (lElmegreen fc Scalol 120041) . The ob served interstellar magnetic field strength is a 



few micro Gauss (e.g. lHeiles fc Trolandll2005l ). in rough equipartition with the kinetic energy 
in the interstellar medium. If the magnetic field were perfectly frozen into the interstellar 
gas during the gravitational collapse of a protostar, the magnetic field strength of a typical 
star like our Sun would be more than 10 orders of magnitude larger than we observe today. 
Thus there must be some mechanis ms that are effectiv e in re moving the excess magnetic flux 



during the star formation process. iMestel &: Spitzerl (119561 ) first suggested that ambipolar 



diffusion (AD) could allow magnetic flux to be redistributed during collapse in low ionization 
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AD-driven gravitational contraction is a quasi-static process. The AD timesca le, £ad, 
in a typical MC is about ten times the free-fall time, tg (IMcKee fc Ostrikerl 120071 ) . Once 
AD has removed a sufficient amount of magnetic flux, a thermally supported core with a 
mass in excess of the Bonnor-Ebert mass will collapse in about a free-fall time. However, 
MCs are observed to be supersonically turbulent. Numerical simulations have shown that 
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turbulence is an efficient mechanism for supporting MCs globally, whi le at the same time 
providing seeds for gravitational collapse by shock compression, (e.g. iKlessen et al.l |2000| ; 



Heitsch et al.ll2001t iLi et al.ll2004l ). By driving large fluctuations in the density, velocity and 



magnetic field , turbulence significantly reduces the AD timescale ( IFatuzzo &: Adams 



2002 



Zweibell 120021 ). AD significantly redistributes magnetic flux. While the importance of this 
process has long been recognized in studies of star formation, it is also important in the 
simpler case in which self-gravity is weak. We therefore wish to determine how AD affects 
the properties of supersonic, magnetized turbulence such as that inferred in MCs. 

It is very challenging to carry out three-dimensional (3D) simulations of AD in molecu- 
lar clouds. The small ionization fraction in molecular clouds means that the ion inertia can 
be neglected. If the AD is treated as diffusion of the magnetic field in a single fluid, the 
time step in explicit codes scales as the square of the grid-size (Ax 2 ), which is prohibitive at 
high resolution; the time step is also proporti onal to the ionization, making it impossible to 
simulate the small ionizations found in MCs (IMac Low et al.lll995l ). Treating the ions and 
neutrals separately as two fluids permits a time step proportional to Ax, but the necessity 
of following the Alfven waves in the ion fluid again leads to very small time steps. Two- 
dimensional fully- implicit codes (e.g. iFiedler fc Mouschoviaslll992l ) were developed to avoid 
this problem, but complex code development would be required to extend this to three di- 
mensions. In addition, implicit treatments can involve multiple iterations to converge, which 
may offset the advantage from the larger timestep. Some attempts have been made to per 
form 3D turbulen ce simulations with AD with semi-implicit schemes (e.g. 
19971 ; lFallell2003l ) but with a heavy cost on computational time. 



Mac Low fc Smith 



To overcome this problem, we introduced the heavy-ion approximation (ILi et al.ll2006l ). 
in which the ionization mass fraction is increased (so as to reduce the ion Alfven velocity) and 
the ion-neutral collisional coupling constant decreased, with the combined result that the 
ion-neutral drag is unchanged. With this approximation, one can perform non-ideal MHD 
turbulence simulations using a two-fluid approach while retaining an accurate treatment of 
the dyna mical interaction be t ween the ions and neutrals in systems with realistic ionization 
fractions. lOishi fc Mac Low! (120061 ) independently made this approximation and used it to 
make a preliminary study of turbulence with AD. Li et al. (2006) discussed the accuracy 
of the heavy-ion approximation and developed criteria for the use of this approximation in 
treating MHD flows with AD. 

In this paper, we investigate the effects of AD on sub-Alfvenic turbulent flows with a 
series of 256 3 and 512 3 MHD turbulence simulations using ZEUS-MPAD (see Li et al. 2006) 
with the heavy-ion approximation. The AD Reynolds number, -Ra o, is an important metric 
in the determination of the significance of AD on turbulent flow ( IZweibel &: Brandenburg 
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19971 ; Li et al. 2006) — for Rad 3> 1, the flow is approximately an ideal MHD flow, whereas 
for i?AD *C 1 the flow of the neutrals approaches a purely hydrodynamic flow. We wish 
to determine the effect of varying the AD Reynolds number on turbulence: How do the 
statistical properties of a turbulent flow depend upon the AD Reynolds number as the 
flow changes from the strong ideal MHD case to the strong AD case? How do the effects of 
turbulent driving of both the neutral and ion gas differ from driving of the neutral component 
alone? What are the criteria necessary to achieve convergence in both the spatial domain 
and in the simulation-time domain when using the heavy-ion approximation? What value 
of the ionization fraction can be used in this approximation and what errors are obtained as 
a function of the ionization fraction used? What is the effect of AD on the velocity and the 
magnetic field power spectra? How do these power spectra compare with recent theoretical 
work on incompressible turbulence in a strong magnetic field? How do the power spectral 
indices compare with more classical turbulent models for smooth, incompressible flows and 
shocked, compressible flows? It is well known that the probability density function (PDF) for 
the density in supersonic isothermal turbulence is log-normal. Is this behavior valid in the 
presence of AD? Finally, how does the presence of ambipolar diffusion in strongly magnetized 
turbulent clouds affect our interpretation of the observed power spectra from these clouds? 

We discuss the heavy-ion approximation and the requirements for its validity in §2. In §3, 
we describe our models based on dimensionless model parameters. Because of the size of the 
parameter space of non-ideal MHD supersonic turbulence, we focus our work on sub-Alfvenic 
turbulence with a thermal Mach number M. = 3. In §4, we report our convergence study 
with the heavy-ion approximation, investigating both spatial and temporal convergence. 
In §5, we report the results on the power spectra of ion and neutral velocities and of the 
magnetic field in turbulence as a function of -Rad- In §6, we discuss the probability density 
function (PDF) of the gas density in the turbulent system. In §7, we investigate other 
physical properties of the turbulent systems that scale with -Rad- We summarize our results 
in §8. Most of the models investigated in details in this paper are based on the assumption 
of ion conservation. However, high density MCs generally have an ionization equilibrium 
timescale that is shorter than the dynamical timescale, so that the ionization inside MCs 
is most likely close to equilibrium. In the Appendix, we demonstrate that nonetheless the 
assumption of ion conservation is generally a satisfactory approximation for molecular clouds. 
The astrophysical implications of our turbulence simulations will be reported in a subsequent 
paper. 
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2. The Heavy-Ion Approximation 

In paper I, we formulated the heavy-ion approximation method for a two-fluid approach 
to MHD simulations with AD. The isothermal MHD equations for the two fluids, ions and 
neutrals, with AD are: 

9pn 



dt 
dpi 

dt 



-V • (p n v n ), (1) 
-V-to-Vi), (2) 



P' n ~Qf = ~Pn(Vn ■ V)V„ - VP„ - 7ADPiPn.(v„ - Vj), (3) 

d v 1 
Pt-^ 1 = -Pi(vi • V)v, - VP; - 7ADPiPn(v, - v n ) + — (VxB) xB, (4) 



dB 

~dt 

V-B = 0, (6) 



Vx(v t xB), (5) 



where p = density, v= velocity, B= magnetic field strength, and 7ad = ion-neutral collisional 
coupling constant (note that in Paper I, this was denoted as 7). The subscripts i and n denote 
ions and neutrals, respectively. Note that in this paper we do not include gravity, so the 
gravitational terms in the two momentum equations have been omitted. In writing these 
equations, we have assumed that ions and neutrals are conserved, as in Paper I. This is 
accurate only if the flow time is small compared to the recombi nation time. Th e inclusion 



of time-dependent chemistry brings in a number of uncertainties ( lDalgarnoll2006l ) ■ However, 
as we show in the Appendix, the contribution from the ionization source terms is in general 
small compared to the AD drag term and can be ignored. 

We define the ionization mass fraction as 

Xi = —, (7) 

Pn 

which we assume to be small. For <C 1, as is the case in molecular clouds, we could 
equally well define x% as the ratio of to the total density p, but in some of our numerical 
models we consider values of x% as large as 0.1, so that p n and p are not equivalent. For 
simplicity, we set the physical value of the ionization mass fraction at a typical observed 
value of x«o,phys = 1CT 6 in our simulations. 



In paper I, we followed iMac Low fc Smith! (119971 ) in implementing a semi-implicit method 



for solving the momentum equations ([3]) and (jlj) in the ZEUS-MP code. The new code, 
ZEUS-MPAD, was tested with several standard AD problems in Paper I. For a two- fluid 
code, the Courant condition restricts the timestep At to be less than Ax/v^i oc y/xi, where 
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v Ai is the ion Alfven velocity. For the low ionizations observed in molecular clouds, Xi ~ 10~ 6 , 
it is still not feasible to perform turbulence simulations even on the latest state-of-the-art 
supercomputing platforms. This problem is compounded because intermittency and the 
highly supersonic nature of the turbulence generate very large contrasts in density. For 
the neutral co mponent, the low density regions could have densities as low as 10~ 4 p n o (e.g. 



Li et al.l 12004 ). and the ion density could be even lower. Therefore, we adopt the heavy-ion 
approximation developed in Paper I, in which the initial ion mass fraction Xio is increased 
but the ion-neutral coupling coefficient 7ad is decreased so as to maintain the same product 
IadXi — 7AD,physXio,ph y s- Following the convention established in Paper I, we set the physical 
value of the ion-neutral coupling coefficient 7AD, P h y s = 9.21 x 10 13 cm 3 g _1 s~ , so that our 
simulations have 7ADXio = 9.21 x 10 7 cm 3 g" 1 s _1 . The heavy-ion approximation reduces 
the frequency of Alfven waves in the ions, which correspondingly increases the integration 
timestep, but it maintains the same dynamical coupling between ions and neutrals. We 
performed three tests in Paper I — the formation of a C-shock, the Wardle instability, and 
a one-dimensional self-gravitating AD collapse — and demonstrated that MHD simulations 
with AD can be sped up by a factor of 10 to 100, depending on the problem, without seriously 
affecting the accuracy. In this paper, we shall consider values of Xio from 10~ 4 to 10 _1 — i.e., 
10 2 — 10 5 times greater than the typical physical value. We shall show that Xio — 10 
corresponding to a speed-up by a factor ~ 100, gives good accuracy. 



-2 



The importance of AD to the flow on a length scale 
diffusion Reynolds number, 



is determined by the ambipolar 



RadW = 



iv 4ir>y AD pip n £v i 



'AD 



where £ad is the AD length scale ( jZweibel fc Brandenburg! 119971 : IZweibell 120021 ) The three 
tests in Paper I all had Rad ~ 1 on the length scale of the problem. Supersonic turbulent 
flows have large contrasts in density, velocity, and magnetic field, and as a result there is a 
large range of length scales involved. The length scales of the local magnetic and velocity 
fields are 



SB 



VSB 

v 



Vv 



(9) 
(10) 



where SB is the change in magnetic field from the mean field B. and where we have as- 
sumed that the mean velocity of the system is zero. Since the ion inertia is negligible in 
the astrophysical problem, it is necessary to ensure that it remains small when the heavy 
ion approximation is used. By comparing the inertia term and the AD drag term in the 
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momentum equations, we deduced in Paper I that this requires 

#AD(4J > ^A» 2 (11) 

where £„ t ~ £sb- We shall verify that this condition is satisfied in our simulations. 

3. Summary of Simulations 

In this paper, we present a series of scale-free turbulence simulations with AD. Three 
dimensionless numbers characterize the simulations: (1) M = 3 1//2 cr nt /c s , the 3D rms Mach 
number of the turbulence, where a nt is the ID rms nonthermal velocity dispersion; (2) the 
plasma (3 = 87rpcj./(B 2 ), which measures the importance of the magnetic field; and (3) 
-Rad(^o) 5 the AD Reynolds number on the scale of the box, which measures the importance 
of AD. The relative importance of the magnetic field on the dynamics of the gas as a whole 
and on the dynamics of the ions is described by the Alfven Mach numbers, 

M A = {(5/2fl 2 M, Mm ~ X ] /2 M A , (12) 

where the expression for Mm is based on the approximation p n ~ p; note that we are 
defining f3 and Ma * n terms of the rms m agnetic field, n ot the mean field. In molecular 



clouds, Ma is observed to be of order unity (ICrutcherlll999l ). whereas Mai is much less than 
unity. 

For -Rad(^o) 3> 1, the dynamics on the scale of the box are described by ideal MHD; if in 
addition, the Alfven Mach number is large (.Ma ^ 1), then the dynamics are approximately 
described by hydrodynamics. In the opposite limit of weak coupling, -Rad(^o) € 1, we expect 
the neutral component to be approximately hydrodynamic, whereas the ions will approximate 
an ideal MHD fluid of their own. Observe that insofar as the neutrals are concerned, the 
cases of low -Rad(^o) and low .Ma could be confused with the case of high -Rad(4>) and high 
.Ma, since in both cases the neutrals behave approximately hydro dynamically. 

The principal goal of this paper is to trace the transition from ideal MHD to weak 
coupling in a turbulent medium by varying the AD Reynolds number of the turbulent box. 
It is not our intention to carry out a complete parameter survey, so we have fixed the plasma- 
(3 = 0.1 and have carried out most of our runs with a thermal Mach number M = 3. The 
corresponding Alfven Mach number is .Ma = 0.67, which is comparable to the observed 
values. We have also carried out a few runs with M = 10, corresponding to A^a = 2.2. 



The simulations are carried out in a cubic box of size £q in each dimension. Periodic 
boundary conditions are applied in all three dimensions, with the intention of approximately 
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representing a small portion of a molecular cloud. The initial magnetic field is oriented along 
the z-axis . We consider the case of turbulence driven according to the recipe described in 



Mac Low! (119991 ): a Gaussian random velocity field with a flat power spectrum in the range 
1 < k < 2, where k = £ /X = /c phys £ /(27r); i.e., the spectrum extends over the wavelength 
range £q > A > £q/2. Random phases and amplitudes are generated in a spherical shell in 
Fourier space and then transformed back into coordinate space to generate each component 
of the driving velocity perturbation. When all three velocity components are obtained, the 
amplitude of the velocity is scaled to a desired initial root-mean-square (rms) velocity, f rms , 
which is defined by the chosen 3D rms Mach number, A4, of the model. Both the ion and 
neutral components start with the same velocity field initially and are driven by a fixed 
driving pattern. We carried out experiments using a variable driving pattern and found 
that the results are statistically indistinguishable from the fixed driving pattern results. We 
also compared driving both the ion and neutral components with driving only the neutral 
component, and again found no statistically significant differences. Therefore, to simplify 
the study, we performed all the simulations using a fixed driving pattern applied to both the 
ions and the neutrals. 

Table 1 lists the initial dimensionless parameters for the models we have calculated. 
The models are labeled as "mxcy", where x = 3 or 10 is the Mach number and y = | log Xio\ 
describes the ionization adopted in the heavy-ion app roximation. For the A4 = 10 models, 



-RadC^o/4) = 1) which is identical to the model used in lOishi &: Mac Low! (120061 ) . We include 
several different values of the AD Reynolds number: .Rad(^o)> which is based on the box size 
and the mean Mach number; (Rad(&v„))vi ^ ne v °l ume average of -Rad based on the neutral 
velocity; and {Rpj)(£ Vi )) v , the volume average based on the ion velocity. For the Mach 3 
models, the latter two agree to within about a factor 2, whereas for the Mach 10 models the 
agreement is within about a factor 3. 



4. Convergence Study 

Simulations of driven turbulence must be converged in both spatial resolution, as is 
the case in any hydrodynamic simulation, and in total simulation time, which is needed to 
reach a steady state. For AD simulations that use the heavy-ion approximation, we must 
also ensure convergence in the mean ionization, xm- As described in §3, we adopt a mean 
physical value for the ionization mass fraction of Xio,ph ys = 10 -6 , but in our simulations we 
use a larger value of Xio and a smaller value of the ion- neutral coupling constant, 7ad, such 
that 7adX«o = 9.21 x 10 7 cm 3 g _1 s _1 is constant. 

The convergence study performed in this section deals solely with globally-integrated 



- 9- 



quantities, such as the total magnetic energy. Figure [T] shows the results of a study of 
Xio- conver g ence ) m which we carried out runs with Xio — 10 _1 , 10~ 2 , 10~ 3 , and 10~ 4 on 
a 256 3 grid with Ai = 3. These four models show a similar evolution pattern, with an 
initial jump in the magnetic energy due to the initial perturbation followed by evolution to a 
quasi-equilibrium with a fluctuating magnetic energy. Fluctuations i n the magnetic energy i n 



AD turbulence have been observed in other simulations as well (e.g. lHawley fc Stondll998l ). 
These fluctuations appear to be random, and they prevent us from carrying out a precise 
convergence study; in particular, we find that runs at different resolutions or with different 
values of Xio yield different time histories of the fluctuations. 

In order to address the issue of convergence in simulation time, we ran the model m3c2 
on a 128 3 grid for a simulation time > 10ij; such a long run is prohibitive using a 256 3 
grid. The flow time is defined as tf = lo/v rms . The time history of the magnetic energy is 
shown in Figure [2j We see that the system is approximately in an equilibrium state after 
one flow time tf, and that the fluctuations persist without a clear period. Since the time 
to reach an equilibrium state is about tf, we continued the Xio-convergence simulations for 
a time somewhat more than 2tf. The mean and standard deviation of the magnetic energy 
determined after the first crossing time are shown in Figure [31 We can see that the total 
magnetic energy converges quickly as Xio decreases. The variation of the mean magnetic 
energy among the models with Xio = 10~ 2 to 10~ 4 is well within the amplitude of the 
fluctuations, but the error in the model with Xio — 10 _1 is larger than this. We conclude 
that the results are converged for Xio ^ 10~ 2 . In Paper I we showed that the heavy-ion 
approximation is satisfied if the AD Reynolds number, -Rad, is large compared to M.m (eq. 
[TTT) . In a turbulent box simulation, it is not trivial to define M.Ai or Rad because of the 
intermittency of turbulence. After driving the box for a period of time, the local A4a{ and 
.Rad have enormous variations — for example, in the Mach 10 model with Xio — 10~ 3 , locally 
defined values of -RadC^J vary by a factor of 10 14 ! Therefore, we compute the volume mean 
(Rad)v f° r both ions and neutrals, using the length scale defined in equation (fTOl) . and list 
them in Table 1. The time at which these volume means are evaluated is also listed in the 
table. Values of (-Rad)v fluctuate, but vary by less than a factor of two for t > tf for both 
the Mach 3 and the Mach 10 models. From Table 1, we see that the requirement to achieve 
a converged solution is actually M-a? / {Rad^v^v ^ 0.03; this is well satisfied for Xio — 
10~ 2 , which accounts for the accuracy of the models with xto — 10 -2 in Figure [31 Note that 
the AD Reynolds number .Rad (A)) evaluated on the scale of the box is significantly larger 
than (-Rad(^J)v, so having A4ai 2 /Rad(^o) "C 1 (or even < 0.03) is not sufficient for the 
validity of the heavy ion approximation. 

To determine the spatial resolution required for turbulent AD simulations, we ran models 
with 128 3 and 512 3 grid cells and with the same initial conditions as model m3c2. The total 
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magnetic energies of these two models are plotted in Figure [3] alongside that of the 256 3 
model. The total magnetic energy, as well as other physical quantities, are converged at a 
resolution of 256 3 using ZEUS-MPAD. 

From these convergence studies, we conclude that we can use the heavy-ion approxima- 
tion with Xio = 10~ 2 to simulate systems with true values of xto ?S 10~ 6 . A spatial resolution 
of at least 256 3 is needed. To obtain reliable statistical results, we suggest driving the system 
for more than 1 one flow time before measuring the physical quantities of the system. 



5. Power Spectra 



The recent work of lOishi fc Mac Low! (120061 ) on MHD turbulence simulations with AD 
compared the magnetic energy spectra with and without ambipolar diffusion and concluded 
that AD produces no dissipation range in the magnetic energy spectrum. In this section, 
we carry out a detailed investigation of velocity and magnetic field power spectra and show 
that AD does in fact have a small, but detectable, effect on the magnetic energy spectrum. 

Generally, for an isotropic turbulent flow, the velocity power spectrum is computed 
using P v (k) = 'Eui(k)u*(k), where Ui(k) is the Fourier transform of the i th component of 
velocity Ui(r) and the sum is over all three velocity components and all wave numbers k 
in the 3D shell k < |k| < k + dk. The inertial range of the power spectrum is expected 
to be a power law P(k) ~ k~ n . For Kolmogorov ( Kolmogorovl 1 1 94 ll ) and Burgers ( Burgers 
19741 ) power spectra, n = 5/3 and 2, respectively. Because of the relatively strong magnetic 
field for {3 = 0.1, especially in our low Mach number models, highly anisotropic distributions 
are expected. It is therefore necessary to compute both P v (k r ) and P v (k z ), the Fourier 
component power spectra perpendicular and parallel to the mean magnetic field, respectively, 
where k r = ^k% + k%. For example, P v (k r ) = Hui(k r )u*(k r ) where the sum is over all three 
velocity components and all wave numbers k r in the 2D shell k r < |k r | < k r + dk r on the 
x — y plane and over all planes along the cylinder with axis z. The magnetic field power 
spectrum is calculated in the same manner. 

Theoretical work on incomp r essible ideal MHD turbulen ce in a strong magnetic field (e.g. 



Goldreich fc Sridharill995l . Il997l : Maron fc Goldreich! l200lh concludes that P v (k 



-5/3 



where k± is perpendicular to the local magnetic field; this has the same exponent as the 
Kolmogorov spectrum. We express this in terms of the exponent in the power spectrum 
as n v i(k±) — 5/3, where we have included t he subscript "i" to indica t e that this applies 



to the ions. However, numerical simulation s (IMaron &: Goldreich! l200ll ; iMiiller et al. 



2003 



Boldyrevll2005l ; iBeresnyak fc Lazarianl 120061 ) give a flatter spectrum that appears consistent 
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with the Iroshnikov-Kraichnan spectrum, n V i(k±) = 3/2 (lIroshnikovlll963l ; lKraichnanlll965l ). 
It is computationally much easier to calculate P v (k r ) than P v (k^), and fortunately the two 
are in close agreement because k r ~ fcj_[l + 0(6 2 )]. On the other hand, P v {k\\) is not 
similar to P v {k z ) be cause k z ~ 6k \, where 6 is the angle between the local magnetic field 



and the ^-direction (IMaron fc Goldreichll200ll ). so we shall not discuss P v {k z ). Since the k 



spectra are often not very meaningful, we also calculate the power spectra in terms of the 
total wavenumber, P{k). These spectra are also of interest because observers do not have 
information on the true direction of the magnetic field. Any measurements of turbulent flows 
inside MCs will be restricted to the line of sight, which can be treated as a random direction 
from the true magnetic field direction. We have demonstrated this by computing the power 
spectra P v (kso) and PB{km) at 60° from the z-axis (the median value of the angle relative 
to the field) for the ideal MHD model m3i, and we find that they agree with the combined 
power spectra P v (k) and Psik) to within the uncertainties. 

In Figure HI we present the time- averaged power spectra of the ion velocity and magnetic 
field for model m3c2 in the time interval ~ (1 — 3)t/. Because of the limited resolution (256 3 ), 
we do not expect the inertial range to extend much beyond k = 10, which corresponds to 
/Cph ys Ax ~ 0.25. Since the driving occurs between k = 1 and 2, we have chosen to infer the 
power-law index by a least-squares fitting of the power spectrum from k — 3 to 10. The 
uncertainty in the index is given by the standard error of the mean, which we calculate as 
the standard deviation evaluated for a total 14 data sets between 1 — 3tf divided by the 
square root of the number of independent samples of the index, which we estimate as 3. 
In order to determine how long it takes for the turbulence system to become uncorrelated, 
we continued models m3i and m3c2 to a time somewhat greater than 5t/. By studying 
the density correlation between data sets at different times, we found that they become 
essentially uncorrelated in a time slightly less than tf. We therefore take the number of 
independent samples to be the largest integer in t Tun /tf. For data sets dumped out between 
1 — 3tf, we shall have 3 independent samples. The range of wavenumbers used to determine 
the power-law index of the power spectrum is very narrow, so we carried out a high-resolution 
run with a resolution of 512 3 (labeled m3c2h) and found that the power-law indexes agreed 
with those from the 256 3 simulation within the errors. We conclude that, although the 
results for the 256 3 runs may not represent accurately the values for the physical case in 
which the inertial range extends over many decades, these results can be used to study the 
dependence of the indexes on the underlying physical parameters. 

No time evolution of the power-law indexes is apparent for t > tf. the time-averaged 
power indexes between tf and 2tf agree with those between 2tf and 3tf within the uncer- 
tainties. Figure UK shows the power spectra, P v ,r{k) and Pe,r(^), of the of the velocity, v r , 
and magnetic field, B r , perpendicular to the global magnetic field. Both ion and neutral 



- 12 - 



velocity spectra are shown. The close agreement between P V i jr {k) and Ps,r{k) is consistent 
with equipartition between th e ion kinetic and magnetic energy perpendicular to the mean 



field (jZweibel &: McKedll995l ; but see §7] below). Figure Sb shows the power spectra of the 
velocity and magnetic field parallel to the global magnetic field, P V)Z (k) and Pe^(fc). The 
power spectra for neutral and ion velocities parallel to the global field are about the same be- 
cause the weakness of the magnetic forces in this direction implies that the ions and neutrals 
are well coupled. Figure shows the combined power spectra P(k) for the neutral and ion 
velocities and for the magnetic field. The power-law indexes resulting from a least-squares 
fit to these spectra are listed in Table 2 and are used to produce the compensated versions 
of these spectra in Figure 0H. 

First, we look at the x«cr conver g ence of the power spectra, using the heavy-ion approx- 
imation. As shown in Table 2, the power-law indexes of the four models m3cl to m3c4, 
which have Xio — 10 _1 — 10~ 4 , are similar within the uncertainties. Figure [5] shows this 
result graphically. Interestingly, the correct value of the power-law index can be obtained 
with Xio = 0.1 even though we found in §4] that convergence in magnetic field energy required 
Xio ^ 0.01. This lack of sensitivity to the value of \%o could be because the relatively low 
resolution of the simulations and the intrinsic fluctuations discussed above lead to significant 
uncertainties in determining the power-law indexes. 

Next, we investigate whether driving the neutrals alone would alter the power spectrum 
of the ion (model m3c2a). In this case, the motion of the ions is mainly due to the drag force 
exerted by the neutrals. The spectral indexes are basically the same as in model m3c2, in 
which both the ions and the neutrals are driven (see Table 2). We conclude that our results 
are insensitive to whether the driving applies to both the neutrals and ions or to the neutrals 
alone. 



How do the spectral indexes depend on the AD Reynolds number? lOishi fc Mac Low 



(120061 ) addressed this question by comparing a run with Rad(^o/^) = 1 to an ideal MHD run, 
both at Mach 10. From visual inspection of their results, they did not find any significant 
difference in the magnetic power spectra [n#(/c) — see their figure 3]. By contrast, when they 
compared a simulation with ohmic dissipation to an ideal MHD simulation, they found a large 
difference in the power spectra. We performed three Mach 10 simulation s with AD (models 



mlOcl to ml0c3) using the same initial conditions as in lOishi fc Mac Low! (120061 ) and a Mach 



10 ideal MHD model mlOi for comparison. The io nization mass fract i on, y^ , in model mlOcl 
is 0.1, which is the same as in the AD model of lOishi fc Mac Low! (120061 ). Unfortunately, 
due to the low densities created in highly supersonic turbulence, it is computationally too 
expensive to continue the Mach 10 models much beyond t = tf. Therefore, we obtained only 
one snapshot of the turbulence for each case of the Mach 10 turbulence with AD, which is 
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insufficient to determine the uncertainty in the indexes. 



We then carried out a quantitative comparison between the AD models and an ideal 
model at Mach 3. The magnetic spectral index for Model m3i, an ideal MHD model with 
the same initial condition as the AD model m3c2, is 71b(&) = 1.25±0.09, which is clearly 
flatter than the value 1.55±0.12 for the AD model. We confirmed this result by comparing 
a high- resolution ideal MHD model (m3ih at 512 3 resolution) with the high- resolution AD 
model m3c2h. However, the effe c t due to AD is much smaller than that of ohmic diffusion 
as reported in lOishi fc Mac Low! (120061 ) . 



More generally, as shown in Table 3, the spectral indexes change systematically with the 
AD Reynolds number. All the indexes listed undergo a statistically significant increase in 
going from the ideal MHD case to the strongest AD case [-Rad (A)) = 0.12]. The indexes for 
the neutral velocity increase down to the lowest value of -Rad, becoming slightly greater than 
2; presumably they level off at yet lower values of R ad, since they approach Burgers value 
of 2 in the hydrodynamic limit (IPadoan et al.l 120071 ). The index for the ^-component of the 
ion velocity [n V i^ z (k)} is locked to that of the neutral velocity since the ions and neutrals are 
well-coupled parallel to the field. With the exception of n vi ^ z , all the ion and magnetic field 
indexes approach constant values at low -Rad, although the value of Rad at which they level 
off varies. At the lowest value of -Rad, the index for the magnetic field has the Iroshnikov- 
Kraichnan value, nn(k) = 1.50 ±0.10, as expected for strong-field (i.e., Iow-A^a) turbulenc e 
jMaron fc Goldreichil200ll ; iMiiller et al]l2003l : iBoldvrevi [2003 : iBeresnvak fc Lazarianl I2OO6} ) . 
However, the effect of the neutrals on the ions is still apparent in this strong AD case, 
since the ion velocity indexes differ significantly from the ideal MHD values. This is to be 
expected, since in the ideal MHD simulation, the power spectrum in the ion fluctuations at 
small scales is due entirely to a cascade from larger scales, whereas in the AD simulations 
the ions are driven at all scales by interactions with the dominant neutrals. Therefore, for 
-Rad(A)) ?S 1, the power spectra of the neutral velocities and of at least the ^-component 
of the ion velocity is close to a Burgers spectrum, and the power spectrum of the magnetic 
field is close to the Iroshnikov-Kraichnan spectrum. 



Passot et al.l (119881) pointed out that the observed linewidth-size scaling, a v oc I 1 / 2 (e.g. 



Solomon et al.l 1 198 71 ). is what would be predicted for Burgers turbulence. IPadoan et al. 



( 200ol ) explained this as the result of ideal MHD turbulence in a weakly magnetized, super- 
sonic (and hence super-Alfvenic) medium. Our results show that the Burgers spectrum can 
be obtained even in sub-Alfvenic turbulence if the AD effect is strong. However, it must 
be borne in mind that our results apply to the inertial range in 256 3 simulations, and as 
discussed above they differ by an unknown amount from the physical case with a far larger 
inertial range. Furthermore, the run with the lowest AD Reynolds number, -Rad(^m) = 0.015, 
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has a value of M.k? / Rad(^vi) comparable to that for the xm = 0.1 runs, which are known 
to be not fully converged; hence, we cannot be sure that this run is fully converged either. 
The trends should be reliable, however. 



6. Probability Density Function 



The probability de nsity function (PDF) for the density of supersonic, isothermal tur- 
bulence is log-normal (jVazquez-Semadenil Il994l ). That is, the volume- weighted or mass- 
weighted probability that the density has a given value is 



f VtM oc exp[-(x ± fi x ) 2 /2al] 



(13) 



where x = ln(p/p), \i x = cr^/2, and the plus a nd minus signs refer to the volume- weighted and 
mass-weighed probabilities, respectively (e.g-. lMcKee fc Ostrikerll2007l ). Hence, the standard 
deviation of the distribution, a x , is related to the means by 

1 



( x )m = 



(14) 



Table 4 gives the values for these quantities and for the median of — x (labeled —x) for a 
range of values of -Rad(^o), extending from ideal MHD [-Rad(-^o) — * °°] to almost decoupled 
[Rad(^o) — 0.12]. All these quantities should be equal for a log-normal PDF, and indeed 
they agree to within the uncertainties, demonstrating that the log-normal behavior of the 
PDF is preserved in the case of ambipolar diffusion. (Note that the equality of the median 
and the mean confirms only that the PDF is symmetric, not that it is a log- normal.) 

Table 4 shows, and Figure M confirms, that the width of the density PDF increases as 
AD becomes more important, although our results do not show a monotonic behavior. This 
increase is plausible due to the decreasing ability of the magnetic field to cushion the shocks 
as .Rad(^o) decreases. For very small -Rad(^o), we expect the dispersion to approach the 
hydrodynamic value, 



a 2 x = ln(l + -M' 



(15) 



(IPadoan fc Nordlundll2002l ). This corresponds to = 0.59 for M. = 3, which is consistent 
with our run at the lowest value of Rad(^o)- F° r the ideal MHD case, IPadoan et al.l (120071 ) 
find a similar r elatio n with the sonic Mach number replaced by the Alfven Mach number. 
Ostriker et al.l (120011 ) do not find such a relation for this case, nor do we. It must be borne 
in mind that the simulations of Ostriker et al. (2001) had a range of values of (3, and that 
our simulations all hav e = 0.1, whi c h is s ubstantially smaller than the value in the super- 
Alfvenic simulations of IPadoan et al.l (120071 ); in addition, our resolution is substantially less. 
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7. Scaling with Rad 

As we have seen in §|5] and $6], the statistical properties of the turbulent box vary as 
the AD Reynolds number of the system, -Rad(^o), changes from the ideal MHD case to the 
strong AD cases. Here we shall examine several other properties as functions of .Rad (A))- 

In the convergence study in §HJ we used the total magnetic energy of the system to gauge 
the convergence in terms of spatial resolution and of ionization mass fraction, xm- Using 
the initial magnetic energy of the system as a reference, we can see how the fluctuating 
magnetic energy of the system, SUb = Ub — Ub , changes as a function of -Rad- Figure [7] 
shows that when Rad is small, 5Ub is also small; indeed, when Rad approaches zero, that 
is when the ion fluid is totally decoupled from the neutral fluid, 5Ub should approach the 
value appropriate for the driven ions (recall that we drive both ions and neutrals). Were we 
to drive only the neutrals, 5Ub would approach zero as Rad — * 0. 

On the other hand, when Rad increases, 5Ub increases to the value appropriate for an 
ideal MHD system. Since MHD waves have equipartition between the kinetic energy normal 
to the field, (1/2) pv±, and the perturbed magnetic energy, SUb, we expect SUb/Ubo = 
(2/3) Ma 2 , under the a ssumption that the velocities are isotropic. Indeed, in their low- (3 



runs, 



Stone et al. ( 19981 ) found SUb/Ubo — 0.6.Ma 2 , consistent with this expectation. For 



our models, A4a = (3M 2 /2 = 9/20, so the theoretical expectation is SUb/Ubo = 0.3. 
We find a significantly smaller value, however, SUb/Ubo — 0.1. We attribute this to our 
boundary conditions: it is possible to have significant kinetic energy that does not perturb 
the field, in the form of edd ies rotating around the field lines or flows along field lines. This 



effect was much smaller for IStone et al.l (119981 ) since they had a much smaller driving scale 



peaked at k — 8. To see whether this effect could be significant in our models, we evaluated 
(v) 2 , where the average is over time and the results are summed over all cells in the box. 
One would expect this to be close to zero, whereas we found it to be a significant fraction of 
(v 2 ). Since these motions are at k ~ 1, however, they do not affect the power spectra. 

Figures [H] and M show density slices for the ions and neutrals normal to the z and y 
axes, respectively, for models m3c2, m3c2rl, m3c2r2, and m3cr3 at t — 3i/. We see that 
the coupling between ions and neutrals gets stronger with larger -Rad- In Model m3c2r3, 
which has -Rad(^o) — 1200 and is very close to ideal MHD, the coupling between the ions 
and neutrals is so strong that there is hardly any difference between the spatial distributions 
of their densities. 

As mentioned above, supersonic turbulence quickly creates large density contrasts in 
the system. In the presence of AD, it also creates large ionization contrasts, as can be 
inferred from Figures [H] and Figures [TUa and [TUb show the ionization mass fraction Xi of 
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a slice at the middle of the turbulence box normal to the z and y axes, respectively, from 
the model m3c2 at time t = 3tf. Regions of low ionization are found to occur in regions of 
high density: Ambipolar diffusion allows shocks to compress the neutrals much more than 
the ions. Because we have assumed overall ion conservation, the mass of ions on a given flux 
tube is constant in the absence of numerical diffusion. Therefore, the change in Xi is purely a 
dynamical result. Furthermore, the contours of logp^ are highly anisotropic and are aligned 
with the z-axis, as shown in Figure [TOb . 

The dispersion in the ionization changes systematically with -Rad- Figure fTTh shows the 
distribution of Xi f° r h ve models, from m3cr-l to m3c2r3, at t — 3i/. We can see thatx; has 
a larger dispersion for smaller values of -Rad- This is to be expected, because in the ideal 
MHD case (-Rad = oo), Xi will maintain a single value for the whole turbulent box, since 
the ions and neutrals are perfectly coupled. With smaller values of Rad, the ions start to 
decouple from the neutrals and the dispersion in Xi increases. We plot the dispersion of x% 
in Figure [Tib. The dispersion decreases as a power of .Rad for models m3c2rl, m3c2r2, and 
m3cr3, which have (-Rad(^J)v > 1, and then becomes approximately constant for the two 
models with (-Rad(^J)v < 1- The turning point is roughly located at (-Rad (A>j ) ) v ~ 1- The 
dispersion in the ionization is somewhat larger than the dispersion in the neutral density, 
even at the smallest values of -Rad we have simulated, since it includes the dispersion in both 
the neutral density and the ion density. 



8. Discussion and Conclusions 

Magnetic fields are an important ingredient in the interstellar medium and are believed 
to play an important role in star formation. On large scales, the magnetic field is frozen 
to the gas, and ideal MHD is appropriate. Indeed, to date almost all 3D simulations of 
MHD turbulence are based on the assumption of ideal MHD. On smaller scales, ambipolar 
diffusion (AD) becomes significant, and in a turbulent medium, the AD lengthscale £ad 
varies substantially due to high contrasts in density, velocity, and magnetic fields. When the 
average value of £ad is comparable to or larger than the size of the turbulent box — i.e., when 
the AD Reynolds number of the box .Rad(^o) ^ 1 — AD can significantly alter the properties 
of the turbulence. Simulating this effect is computationally challenging, however, since the 
timestep required in explicit codes is proportional to both the square of the gridsize, Ax 2 , 
and to the square root of the ionization mass fraction, Xio, both of which are exceedingly 
small for accurate modeling of MCs. Semi-implicit treatments can avoid the problem with 
Ax 2 but not the one with the small ionization mass fraction. To overcome the latter problem, 



Oishi fc Mac Low! (120061 ) adopted an artificially high value for the ionization mass fraction, 
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Xm = 0.1, and a correspondingly low value for the ion-neutral coupling coefficient, 7 ad, in 
carryi ng out the 3D simulation for a turbulent medium with ambipolar diffusion. iLi et al. 



( 120061 ) independently developed this approximation, which they termed the heavy-ion ap- 
proximation, and determined the condition for its validity by testing it on several classical 
MHD problems involving AD. 

In this paper, we report the results of our simulations of sub-Alfvenic turbulence with 
AD using the heavy-ion approximation. We assume that the ions are conserved, but in the 
Appendix we show that our results also apply approximately to the case of time-dependent 
ionization for realistic molecular densities. Our models focus on the case of a thermal Mach 
number of 3 and a plasma (3 of 0.1, corresponding to an Alfven Mach number M.a = 0.67. 
By using this relatively low value of the Mach number, we are able to perform a number of 
256 3 simulations with a duration of 3— 5tf, where tf = £o/v rm s is th e flow time across the box. 



We find, in agreement with previous workers (e.g. lHawley fc Stondll998l ). that simulations 
of turbulence with AD have significant fluctuations in most physical quantities, which makes 
it difficult to accurately determine the statistical properties of the system. We carried out 
several convergence studies to determine the validity of our simulations: First, we showed 
that independent samples of the turbulence can be obtained at time intervals — tf beginning 
at t — tf, and that a total running time of ~ 3tf is sufficient. Second, we showed that the 
results are converged to within the uncertainties for a spatial resolution of 256 3 . Finally, we 
showed that the heavy-ion approximation with x%o = 10 -2 can represent a xm — 10~ 6 system 
as observed in MCs with sufficient accuracy. This speeds up the calculation by a factor of 
about 100, but it is nonetheless a factor 10 slower than an ideal MHD simulation. 

High-resolution ideal MHD turbulence simulations show that the velocity power spec- 
trum for super-Alfvenic turbulence will be close to a Burgers spectrum with a power index 
n v (k) ~ 2 (Padoan et al. 2007). If the turbulence is sub-Alfvenic, the power spectrum will 
be highly anisotropic with respect to the direction of the magnetic field. Theoretical studies 
on incompressible turbulence suggest the velocity power spectrum no rmal to the magnetic 
field will be a Kolmogorov-like spectrum, with power index of 5/3 (e.g. iGoldreich &: Sridhar 
19951 ). Numerical simulatio ns indicate that in strong fields the power in dex is close to 1.5 
(i.e., at low values of M A ) JMaron fc Goldreichl lioOll : iMuller et al.ll2003h . and that the 5/3 
index is realized only for relatively weak mean fields ( Boldyrev 12005 ). 



We have studied how the power spectra change as a function of the importance of 
ambipolar diffusion, which is measured by the AD Reynolds number -Rad- We use two 
different values of the AD Reynolds number: -Rad (A)) is defined in terms of the box size and 
the rms velocity in the box, whereas (Rab(K 1 ))v is the volume average of the AD Reynolds 
number defined in terms of local parameters, with a length equal to the scale over which 
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the ion velocity varies. (-Rad(^J)v is the quantity that enters the criterion to ensure the 
validity of the heavy- ion approximation (eq. ITTT) . It varies with time during a simulation, 
but is constant to within a factor of 2 in the simulations reported here; once the system 
reaches equilibrium at t ~ tf, the change in (-Rad)v is < 10%. For the five models with 
(-Rad(^J)v from 0.015 to 215.83 [initial .Rad (A)) from 0.12 to 1200] that we have computed, 
we see a progressive transition of system properties from a model with a strong AD effect 
to a near ideal MHD model; we also computed ideal MHD models for comparison. All the 
power law indexes we computed increase in going from the ideal MHD case to the case of 
strongest AD, and most of them appear to approach a constant at small -Rad- 

For the ideal MHD case, we confirm that the power spectrum of the ion velocity normal 
to the field is consistent with the Iroshnikov-Kraichnan spectrum (n vi ^ r = 1.5), as found in 
previous studies of turbulence in strong fields. This index increases with the importance of 
AD; it is consistent with the Kolmogorov value in four of our AD models, but is larger for 
the case in which the AD is strongest. The perpendicular magnetic field power is usually 
larger than that of the parallel magnetic field power, and as a result the power-law index 
for the total field, riB(k), is about the same as that for the perpendicular components of 
the field, riB >r {k). We find that these power-law indexes are about 1.2 in the case of ideal 
MHD and rise to about 1.5 for the case in which AD is strongest. This does not agree with 



the conclusion of lOishi &: Mac Low! (120061 ). who found no difference in the magnetic power 
spectra of an ideal MHD model and models with strong AD. We note that the change in the 
index of the power spectra we find between ideal and AD-dominated MHD is small compared 
to the difference they reported between ideal MHD and MHD dominated by ohmic diffusion. 

By comparing the volume- averaged value of In p, the mass- averaged value, the dispersion 
in values, and the median, all of which have equal magnitude for a log-normal PDF, we 
concluded that the density PDF is indeed log-normal for all the cases we considered. The 
dispersion increases system atically as AD increases in importance, and is consistent with the 



Padoan fc Nordlundl (120021 ) result for the strong AD cases. 



An important result from our sub-Alfvenic turbulence simulations with AD is that the 
neutral gas in systems with small f3 (strong magnetic field) and strong AD (small -Rad) 
behaves like that in systems with large f3 (weak magnetic field) and no AD. In particular, 
the neutral-velocity power spectrum in a strongly magnetized medium with strong AD is 
approximately consistent with a Burgers spectrum [n v (k) ~ 2]. It is thus not possible to 
infer the strength of the magnetic field from observations of the power spectrum unless it is 
known that the observations are on a sufficiently large scale that AD is not important. 
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A. ION CONSERVATION VERSUS TIME-DEPENDENT IONIZATION 



The calculations we have discussed in this paper are based on the assumption that the 
number of ions is conserved. In fact, as the density changes, ionization and recombination 
will change the number of ions. The ionization timescale is 

'3 x KT 17, 



100 



e vm-7;v - I yr ' (A1) 

S>CR VJ -U / \ (,cr / 

where x e = n e /% is the ionization number fraction and Ccr is the i onization rate per H 
atom, which is inferred to be about (2.5 — 5) x 10~ 17 s _1 in dense clouds ( lDalgarnoll2006l ) . The 



recombination time scale is t Tec = l/an e , where a is the relevant recombination coefficient. 
In equilibrium these two time scales are equal, which implies that the equilibrium ionization 
is 
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where the numerical eva luation is for £cr = 3 x 10~ 17 and a = 2.5 x 10 
( IMcKee &: Ostrikerl 120071 ). [This estimate of the ionization is based on the assumption that 
HCO + dominates the ionization; if small PAHs dominate, then the effective recombination 



rate is about 10 times smaller (IWakelam fc Herbstl 120081 ) and the equilibrium ionization is 
several times larger.] Furthermore, one can show that if the ionization is close to equilibrium, 
then the the e-folding time for the ionization to approach equilibrium is half as large as the 
ionization and recombination times: 
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Note that the ionization time scale is generally ord ers of magnitude less than the chemical 
equilibration time scale, which can be ~ 10 5 yr (e.g. iPadoan et al.ll2004J : IWakelam fc Herbst 
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20081 ). In a molecular cloud with low ionization, the ionization time scale is short compared 



to the typical dynamical time scale, 

R R 



^dyn 



o 0.72RH 2 km s- 1 PC 



1.36 x 10 6 i?^ 2 yr, (A4) 



where R pc = R/(l pc) and whe re we have assumed that the velocity dispersion obeys the 
standard linewidth-size relation (IMcKee fc Ostrikerl 120071 ). Gas in molecular clouds is there- 
fore expected to be close to ionization equilibrium except in regions where the dynamical 
time scale is short, as in shocks. 

One can include ionization and recombination as source terms in the mass and momen- 
tum equations (jTJ) — (TJJ as: 

>nV„) - Si, (A5) 

HvJ + Sx, (A6) 

(A7) 

(VxB)xB + 5 2 . (A8) 
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•Si = (CcrUh - an 2 ) rrn , (A9) 



an 2 



S 2 = Cc^uum |^v n - - — ^-Vij , (A10) 

where rrii is the ion mass and we have assumed charge neutrality, rii = n e . Since the 
ionization generally low, much of the ionization of H2 will be transferred to heavy molecules 
such as HCO + . Therefore, we use m.j also i n the ionization compon e nt of the source terms 
Si and S2; this differs from the treatment in iBrandenburg fc Zweibell (119951 ). In equilibrium, 
CcR,nf/ = «^e, eq ; so ^ ne coefficient of Vj in equation (lAlOj) is simply (n e /n e ^ eq ) 2 . As a result, 
the second term in the momentum source term dominates when the gas is overionized. 

Let TZi on be the ratio of the momentum source term S2 to the AD drag term. Using 
equation flAll) . we find that in equilibrium this ratio is 

T^-ion.cq = = " T — • (AH) 

We see that the ionization/recombination source term is important only when the density is 
very low and/or the ionization timescale is very small. However, in MCs these conditions are 
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generally not satisfied, so that TZi on is very small and the ionization/recombination source 
terms can be ignored. For example, the t ypical density and ioniza tion in MCs are nn,3 = 
?t,h/(10 3 cm -3 ) > 1 and x e ~ lO -7 ^ 1 ^ 2 ( McKee fc Ostriker 2007 ). so that the ionization 
timescale is t eq ~ 60n H 1 g 2 yr and 7^i on ,eq ~ 1.3 x 10~ 3 n H 3 . 

To verify that the momentum source terms are indeed negligible in realistic cases, we 
modified ZEUS-MPAD to include the source terms and performed a model simulation based 
on the initial conditions of model m3c2 (which has 7^ ioni e q = 0) but with 7^ ion eq = 1.7 x 10 -3 . 
This corresponds to a density nn = 550 cm -3 ; since this is smaller than the typical density in 
molecular gas, this represents an approximate upper bound on the effect of time-dependent 
ionization. The mean value of 7^i onjeq in the model is about a factor of 1.6 times the initial 
value using equation fl Al 1 [) : 7£ ion>cq oc 1 /pi/ 2 , and for a log-normal distribution one can 
show that the mean of (pn/ Pn) 1 ^ 2 is exp(3cr 2 /8), where a 2 is the dispersion of the log normal. 



Padoan fc Nordlundl (J2002J) estimate a 2 ~ ln[l + (M/2) 2 }, which is 1.18 for our Mach 3 



models; hence, (7£ ioilieq ) ~ 2.8 x 10~ 3 . However, non-equilibrium effects are very important: 
For those cells that are overionized, TZ lon will be larger by a factor of order an 2 / CcR n H = 
(n e /n eeq ) 2 . The average value of {n e /n e!eq } 2 is about 16, and as a result the average value of 
7li on is 0.02, significantly larger than the equilibrium value. Although the mean value of 7Zi on 
is small, time-dependent ionization has a detectable effect on the spectra of the turbulence, 
being about 1 a steeper than those for the case of ion conservation. For realistic molecular 
densities, 7^ ion will be smaller, so the spectra for the time-dependent case will be closer to 
those for the conservation case. We conclude that MC models with time-dependent ionization 
are generally well approximated by models using the assumption of ion conservation. 



REFERENCES 

Beresnyak, A., & Lazarian, A. 2006, ApJ, 640, L175 
Boldyrev, S. 2005, ApJ, 626, L37 

Brandenburg, A., & Zweibel, E. G. 1995, ApJ, 448, 734 

Burgers, J. M. 1974, The Nonlinear Diffusion Equation (Dordrecht: Reidel) 

Crutcher, R.M. 1999, ApJ, 520, 706 

Dalgarno, A. 2006, PNAS, 103, 411 

Elmegreen, B. G. & Scalo, J. 2004, ARA&A, 42, 211 

Falle, S. A. E. G. 2003, MNRAS, 344, 1210 



-22 - 



Fatuzzo, M. & Adams, F.C. 2002, ApJ, 570. 210 

Fiedler, R. A. & Mouschovias, T. Ch. 1992, ApJ, 391, 199 

Fiedler, R. A. & Mouschovias, T. Ch. 1993, ApJ, 415, 680 

Goldreich, P. & Sridhar, H. 1995, ApJ, 438, 763 

Goldreich, P. & Sridhar, S. 1997, ApJ, 485, 680 

Hawley, J. F. & Stone, J. M. 1998, ApJ, 501, 758 

Heiles, C. & Troland, T. H. 2005, ApJ, 624, 773 

Heitsch, F., Mac Low, M.-M., & Klessen, R. S. 2001, ApJ, 547, 280 

Iroshnikov, P. S. 1963, AZh, 40, 742 (English transl. Soviet Astron., 7, 566 [1964]) 

Klessen, R. S., Heitsch, F., & Mac Low, M.-M. 2000, ApJ, 535, 887 

Kolmogorov, A. 1941, Dokl. Akad. Nauk SSSR, 31, 538 

Kraichnan, R. H. 1965, Phys. Fluids, 8, 1385 

Li, P. S., Norman, M. L., Mac Low, M.-M., & Heitsch, F. 2004, ApJ, 605, 818 
Li, P. S., McKee, C. F., & Klein, R. I. 2006, ApJ, 653, 1280 
Lizano, S. & Shu, F. H. (1989), ApJ, 342, 834 

Mac Low, M.-M., Norman, M. L., Konigl A., & Wardle, M. 1995, ApJ, 442, 726 

Mac Low, M.-M. & Smith, M. D. 1997, ApJ, 491, 596 

Mac Low, M.-M. 1999, 524, 169 

Maron, J. & Goldreich, P. 2001, ApJ, 554, 1175 

McKee, C. F., & Ostriker, E. C. 2007, ARAA, in press. 

Mestel, L., & Spitzer, L. 1956, MNRAS, 116, 503 

Mouschovias, T. Ch. 1976, ApJ, 207, 141 

Mouschovias, T. Ch. 1977, ApJ, 211, 147 

Mouschovias, T. Ch. 1979, ApJ, 228, 475 



-23- 

Miiller, W. C, Biskamp, D., & Grappin, R. 2003, Phys. Rev. E, 67, 066302 

Nakano, T., & Nakamura, T. 1978, PAS J, 30, 671 

Nakano, T. & Tademaru, E. 1972, ApJ, 173, 87 

Oishi, J. S. & Mac Low, M.-M. 2006, ApJ, 638, 281 

Ostriker, E. C, Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980 

Padoan, P. & Nordlund, A. 2002, ApJ, 576, 870 

Padoan, P., Zweibel, E., & Nordlund, A. 2000, ApJ, 540, 332 

Padoan, P., Willacy, K., Langer, W. & Juvela, M. 2004, ApJ, 614, 203 

Padoan, P., Nordlund, A., Kritsuk, A. G., Norman, M. L., & Li, P. S. 2007, ApJ, 661, 972 

Passot, T., Pouquet, A., & Woodward, P. 1988, A&A, 197, 228 

Shu, F. H. 1983, ApJ, 273, 202 

Solomon, P. M., Rivolo, A. R., Barret, J. & Yahil, A. 1987, ApJ, 319, 730 
Spitzer, L., Jr. 1968, Diffuse Matter in Space (New York: Interscience) 
Stone, J.M., Ostriker, E. C, & Gammie, C. F. 1998, ApJ, 508, L99 
Vazquez-Semadeni, E. 1994, ApJ, 423, 681 

Wakelam, V., & Herbst, E. 2008, ArXiv e-prints, 802, arXiv: 0802 .37571 

Zweibel, E. G. 2002, ApJ, 567, 962 

Zweibel, E. G. & Brandenburg, A. 1997, ApJ, 478, 563 

Zweibel, E. G., & McKee, C. F. 1995, ApJ, 439, 779 



This preprint was prepared with the AAS IATgX macros v5.2. 



-24- 




Fig. 1. — Time evolution of the total magnetic energy, Ub, normalized to the initial total 
magnetic energy Ub,o, for models m3cl (xm — 10 _1 , solid line), m3c2 (xm = 10~ 2 , dash line), 
m3c3 (xm = 10~ 3 , dot-dash line), and m3c4 (xw — dotted line). The systems settle 

into approximate equilibrium states for t > tf. 



-25 - 




Fig. 2. — Time evolution of the total magnetic energy, Ub, for four 128 3 turbulence models 
with the same initial conditions (models m3cl to m3c4). One of the models (m3c2, with 
Xio = 1CT 2 ) runs until t > 10t/. The system is approximately in equilibrium for t >tf, with 
random fluctuations. 
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Fig. 3. — Convergence behavior of the time-averaged total magnetic energy, (Ub), for models 
m3cl to m3c4 as a function of Xio- The total magnetic energies (circles) are averaged after 
the first crossing time and the error bars show the standard errors of the means. The total 
magnetic energy is converged within the fluctuation limits for an ionization mass faction 
Xio < 1CT 2 . The total magnetic energy of two models with the same initial conditions of the 
model m3c2 but with resolution of 128 3 (diamond) and 512 3 (square) are also plotted. The 
Xm of these two models are the same as m3c2 but changed here by a small amount in the 
plotting for the clarity of the overlapping error bars. 
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Fig. 4. — (a) Velocity power spectra of the neutrals (n TOjT ., solid line) and ions (n vijr , dashed 
line) and the magnetic field power spectrum (ng !r , dot-dash line), all perpendicular to the 
global magnetic field for model m3c2. (b) Same as (a) but the components are parallel 
to the global magnetic field. The ions and neutrals have very similar spectra parallel to 
the global magnetic field direction since only weak fields are induced in the perpendicular 
direction, (c) The combined 3D velocity power spectra of neutrals (n vn , solid line) and ions 
(n V i, dashed line), and the power spectrum of the magnetic field (ng, dot-dash line), (d) 
The compensated 3D velocity power spectra of neutrals (solid line) and ions (dashed line), 
and the compensated magnetic field power spectrum (dot-dash line). The power law indexes 
used for the compensation are listed in Table 2. 
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Fig. 5. — The compensated power spectra of models m3cl (xm = 10 _1 , solid line), m3c2 
(Xio = 10~ 2 , dashed line), and m3c3 (xio — 10~ 3 , dot-dash line). The spectra are compen- 
sated by the power law indexes of the inertial range fitted between k — 3 — 10. The thin 
lines are neutral velocity spectra and the thick lines are the magnetic field power spectra. 
The spectra are shifted up and down for side-by-side comparison. As shown in the figure 
and listed in Table 2, the power law indexes are not sensitive to the choice of ionization mass 
fraction Xio, even for Xio as large as 0.1. 
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Fig. 6. — Density PDF of ideal MHD model m3i (circles) and model m3c2 (squares). The 
density PDF of model m3c2 shows significantly larger dispersion than the ideal MHD model 
because of the effects of ambipolar diffusion. The dispersion and mean of the PDFs are listed 
in Table 4. 
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Fig. 7. — Time- averaged change in the normalized total magnetic energy, SUb/Ubo for five 
models m3c2r-l to m3c2r3 as a function of Rad^vJ- Uncertainties are shown as error bars. 
The dashed line is 8U B /U B o for the ideal MHD model m3i. With increasing Rad, SU b /U B o 
approaches the ideal MHD model value. 
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Fig. 8. — Logarithmic density (logp) slices of models m3c2 (1st row), m3c2rl (2nd row), 
m3c2r2 (3rd row), and m3c2r3 (4th row) at the middle of the turbulent box normal to the 
z-direction at time t = 3tf. The left column shows the neutral density and the right column 
shows the ion density. When (Rad)v is large, the ions and neutrals are sufficiently strongly 
coupled that they evolve like a single fluid. 




Fig. 9. — Same as Figure [H] but the slices are at the middle of the turbulent box normal to 
the ^-direction. The ion density is highly anisotropic due to the strong magnetic field and 
relatively weak turbulence (.Ma < 1)- 
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Fig. 10. — Spatial distribution of ionization mass fraction Xi- ( a ) The Xi °f a snce &t the 
middle of the turbulence box normal to the z-direction (B- field direction). The contours are 
log Xi an d the grey scale (color scale in online version) map is log p n . Small Xi regions usally 
associate with high density regions, (b) Same as (a) but the slice is normal to the y-direction. 
The contours are highly anisotropic because of the restraint of ions due to strong magnetic 
field. 
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Fig. 11. — i?AD effects: (a) The PDFs of \i f°r models m3c2r-l (solid), m3c2 (dash), 
m3c2rl (dotted), m3c2r2 (dot-dash), and m3c2r3 (thick dash) at t = 3tf as a function 
of (-Rad(AjJ)v- (b) The time-averaged dispersions of the x% distributions for the five models 
versus (-Rad(4>J)v ov er 2 tf. The x% an d (-Rad(4>J)v show a power law relation when the 
(-Rad(^J)v > 1- When (-Rad(^J)v < 1, the dispersion of x% approaches a constant. 
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Table 1: Model Parameters and Statistical Results 



Model* 


M rms 


Xi 


7AD 


Time 

(*/) 


-Rad(^o) 






M 2 Ai /(R A n(£ v M 
(xl0~ 3 ) 


m3cl 


3 


io- 1 


4 


2 


1.2 


0.3916 


0.2184 


167.5 


m3c2 


3 


io- 2 


40 


3 


1.2 


0.3206 


0.2197 


16.01 


m3c3 


3 


io- 3 


400 


2 


1.2 


0.3392 


0.1908 


1.849 


m3c4 


3 


10~ 4 


4000 


2 


1.2 


0.3363 


0.1894 


0.213 


m3c2a** 


3 


io- 2 


40 


3 


1.2 


0.3374 


0.2183 


15.17 


m3c2h** s 


3 


io- 2 


40 


3 


1.2 


0.2709 


0.1750 


20.17 


m3c2r-l 


3 


io- 2 


4 


3 


0.12 


0.0298 


0.0145 


148.1 


m3c2rl 


3 


io- 2 


400 


3 


12 


2.9266 


2.5480 


1.767 


m3c2r2 


3 


io- 2 


4000 


3 


120 


25.455 


25.090 


0.179 


m3c2r3 


3 


io- 2 


40000 


3 


1200 


227.47 


225.59 


0.020 


mlOcl 


10 


io- 1 


4 


1.25 


4 


1.0720 


0.4326 


909.6 


ml0c2 


10 


IO" 3 


40 


1.25 


4 


1.1596 


0.4172 


90.24 


ml0c3 


10 


IO" 4 


400 


1.25 


4 


1.2653 


0.4249 


8.860 


m3i^ 


3 


oo 


oo 


3 


oo 


oo 


oo 




m3iht*** 


3 


oo 


oo 


3 


oo 


oo 


oo 




mlOit 


3 


oo 


oo 


3 


oo 


oo 


oo 




* Models 


are labeled as 1 


: 'mxcy," 


where x 


is the thermal Mach number and y = 


1 logXiol- 



Models labeled "mxcyrn" have i?AD(4) = 1.2 x 10 n . 

** Driving applied only to the neutrals. 

*** High resolution model (512 3 ). 

t Ideal MHD models. 

^ Root mean squared (rms) values. 



Table 2. Spectral indexes of Velocity and Magnetic Field Power Spectra for Models in Convergence Studies 







n ■ (k) 






n (k) 
' t vn,zV /i I 


n (k\ 
< i vn \">) 


r>r? (k) 






m3cl 


1 


87+0 1Q 

• Kj t —1—\J . ± U 


1.87±0.09 


1.87±0.10 


1.99±0.06 


1.99±0.10 


1.99±0.06 


1.41±0.16 


2, 


98+0 iq 


1.55±0.16 


m3c2 


1 


79+0 14 

. 1 <J —1—\J . ± t: 


1.92±0.07 


1.88±0.08 


1.92±0.07 


1.90±0.07 


1.91±0.06 


1.44±0.12 


2, 


.12+0.13 


1.55±0.12 


m3c3 


1 


,75±0.18 


1.91±0.09 


1.86±0.07 


1.94±0.04 


1.90±0.09 


1.93±0.04 


1.40±0.12 


2. 


,02±0.18 


1.50±0.13 


m3c4 


1 


.79±0.16 


1.90±0.06 


1.82±0.11 


1.91±0.07 


1.90±0.07 


1.91±0.05 


1.42±0.12 


2 


.08±0.16 


1.53±0.13 


m3c2h 


1 


,64±0.06 


1.98±0.03 


1.86±0.03 


1.92±0.03 


1.95±0.03 


1.94±0.03 


1.48±0.05 


2. 


.14+0.07 


1.56±0.05 


m3i 


1 


,46±0.10 


1.31±0.10 


1.41±0.08 








1.17±0.09 


1 


,72±0.11 


1.25±0.09 


m3ih 


1 


,48±0.06 


1.38±0.05 


1.45±0.05 








1.14±0.07 


1. 


,61±0.08 


1.23±0.07 


m3c2a 


1 


,83±0.14 


1.93±0.07 


1.90±0.06 


1.93±0.05 


1.90±0.07 


1.92±0.04 


1.40±0.10 


2 


,03±0.16 


1.50±0.10 


mlOcl 




1.50 


1.73 


1.57 


2.00 


1.79 


1.94 


0.99 




1.29 


1.05 


ml0c2 




1.05 


1.86 


1.34 


2.00 


1.93 


1.98 


1.11 




1.63 


1.21 


ml0c3 




1.05 


1.73 


1.34 


2.04 


1.84 


1.98 


1.15 




1.19 


1.16 


mlOi 




1.11 


1.53 


1.27 








1.38 




1.18 


1.34 
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Table 3: Spectral indexes of Velocity and Magnetic Field Power Spectra for M = 3 Models 
with Varying Rad 





m3c2r-l 




m3c2 




m3c2rl 


m3c2r2 


m3c2r3 




m3i 


-Rad(^o) 


0.12 




1.2 






12.0 




120 


1200 




oo 


(R A D{L t )) V 


0.015±0.002 


0. 


,21±0. 


.01 


2 


,74±0. 


.13 


25.25±0.79 


215.83±5.00 




oo 




1.89±0.15 


1 


,79±0. 


14 


1. 


,78±0. 


13 


1.72±0.12 


1.67±0.10 


1 


,46±0.10 




2.07±0.06 


1 


,92±0. 


.07 


1. 


,80±0. 


.14 


1.32±0.12 


1.27±0.16 


1 


,31±0.10 


n vi (k) 


1.92±0.08 


1 


,88±0. 


.08 


1 


,82±0. 


.12 


1.56±0.08 


1.50±0.09 


1 


,41±0.08 




2.17±0.05 


1 


,92±0. 


.07 


1 


,94±0. 


.09 


1.85±0.12 


1.67±0.10 








2.07±0.06 


1 


.90±0. 


.07 


1 


.80±0. 


.07 


1.32±0.12 


1.27±0.16 






n VD (k) 


2.14±0.04 


1 


,91±0. 


06 


1 


.87±0. 


.09 


1.57±0.15 


1.50±0.09 






n B ,r( k ) 


1.45±0.10 


1 


,44±0. 


.12 


1 


.38±0. 


.15 


1.27±0.13 


1.12±0.08 


1. 


.17±0.09 


n B , z (k) 


2.03±0.12 


2 


,12±0. 


13 


2 


,08±0. 


.12 


2.03±0.13 


1.82±0.08 


1. 


,72±0.11 


n B (k) 


1.53±0.09 


1 


,55±0. 


.12 


1. 


,50±0. 


.14 


1.38±0.12 


1.22±0.08 


1 


,25±0.09 


n vi(k r ) 


1.95±0.09 


1 


.99±0^ 


06 


1 


.96±0. 


13 


1.64±0.08 


1.54±0.09 


1. 


.60±0.07 


^vn ( k r ) 


2.11±0.04 


2 


.08±0. 


.06 


1 


.83±0. 


.09 


1.66±0.08 


1.55±0.09 






n B (k r ) 


1.50±0.10 


1 


,53±0. 


.11 


1. 


,36±0. 


.15 


1.29±0.13 


1.29±0.08 


1 


,33±0.09 
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Table 4: Statistical Parameters of the Density PDF 



Model 


(RAB(£ Vi )) V 


~( x )v 


( x ) M 


—x 


¥1 


m3c2r-l 


0.015±0.002 


0.58±0.03 


0.55±0.02 


0.56±0.04 


0.60±0.03 


m3c2 


0.21±0.01 


0.65±0.04 


0.60±0.04 


0.62±0.05 


0.70±0.05 


m3c2rl 


2.74±0.13 


0.59±0.07 


0.51±0.05 


0.49±0.07 


0.70±0.12 


m3c2r2 


25.25±0.79 


0.32±0.03 


0.32±0.03 


0.33±0.03 


0.32±0.02 


m3c2r3 


215.83±5.00 


0.39±0.03 


0.37±0.03 


0.38±0.05 


0.40±0.02 


m3c2h 


0.18±0.01 


0.60±0.03 


0.56±0.02 


0.58±0.02 


0.64±0.05 


m3i 


oo 


0.43±0.04 


0.40±0.04 


0.40±0.07 


0.46±0.04 


m3ih 


oo 


0.40±0.05 


0.38±0.06 


0.39±0.09 


0.42±0.05 



